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We consider a dilute gas of hard spheres under shear. We use one of the predominant models to 
study this system, namely the SLLOD equations of motion, with an iso-kinetic Gaussian thermostat 
in between collisions, to get a stationary total peculiar kinetic energy. Based on the previously 
obtained result that in the non-equilibrium steady state and in case the number of particles -/V 
becomes large, the coefficient of dynamical friction representing the iso-kinetic Gaussian thermostat 
for the SLLOD dynamics fluctuates with l/ViV fluctuations around a fixed value, we show on 
analytical grounds that for a hard sphere gas at small shear rate and with a large number of 
spheres, the conjugate pairing of the Lyapunov exponents is expected to be violated at the fourth 
power of the constant shear rate in the bulk. 



PACS Numbers: 05.20.-y, 05.45.-a, 05.60. Cd 



I. INTRODUCTION 



Non-equilibrium molecular dynamics (NEMD) simula- 
tions of Navier-Stokes equations have been used to study 
the shear viscosity properties of fluids for a long time. To 
study the coefficient of shear viscosity in Navier-Stokes 
equations, a carefully chosen periodic boundary condi- 
tion in NEMD simulations is enough to drive the system 
out of equilibrium. Based on these ideas, in the beginning 
days of the development of this subject, an algorithm was 
constructed from simple Newtonian equations of motion 
using a periodic boundary conditions, the so-called Lees- 
Edwards boundary conditions 0. However, it was soon 
realized that in the absence of an explicit dependence on 
the shear field in this algorithm, one could not make an 
appropriate connection with the Green-Kubo relations, 
and therefore, it was difficult to deal with from a theo- 
retical point of view. As a remedy, some other algorithm 
with an explicit dependence on the shear field was called 
for, and the DOLLS and the SLLOD algorithms were 
born. 

The main idea behind the DOLLS and the SLLOD al- 
gorithm is an explicit dependence on the shear field, 7. 
The DOLLS algorithm was implemented first The 
SLLOD equations of motion were proposed soon after 
||, and are now preferred because they are equivalent 
to the boundary driven method (J]. Both algorithms 
have to be supplemented by a thermostat, which con- 
tinuously removes the energy generated due to the work 
done on the system by the shear field such that an non- 
equilibrium steady state (NESS), homogeneous in space, 
can be reached. 

In this paper, we will look at the SLLOD equations 
of motion for a gas of hard spheres from the point of 



view of dynamical systems. The Lyapunov exponents 
of a system of particles, obeying the SLLOD equations 
of motion and mutually interacting by means of WCA 
potential with an iso-kinetic Gaussian thermostat, was 
first studied by Morriss ||. The study showed that the 
shear viscosity can be obtained from the sum of all the 
Lyapunov exponents. The simulation results in Ref. j^] 
also indicated that once the Lyapunov exponents are ar- 
ranged in ascending order of magnitude, the sum of the 
largest and the smallest, the sum of the second largest 
and the second smallest and so on, were the same. The 
phenomenon of such pairing of the Lyapunov exponents 
is known as the Conjugate Pairing Rule, or the CPR. 
Since it is in general difficult to calculate all the individ- 
ual Lyapunov exponents of a system, an extensive the- 
oretical study soon ensued to understand the CPR for 
the Lyapunov exponents of systems obeying the SLLOD 
equations of motion. Evans and co-workers investigated 
this point Q to conclude that the Lyapunov exponents 
pair exactly for general inter-particle potentials and all 
7. In a follow-up work, Sarman et al. carried out 
simulation studies in support of Ref. H . 

In the next few years, the connection between the dy- 
namical systems theory and statistical mechanics saw a 
surge of interest. Some situations were found, where it 
waspossible to prove that the CPR is satisfied exactly 
§-|fl). The status of the CPR for the SLLOD and the 
DOLLS dynamics was revisited. For a system of parti- 
cles obeying the SLLOD and the DOLLS dynamics with 
a WCA inter-particle potential and arbitrary 7, CPR was 
reported to be violated on the basis of simulation results 
fll2| , [i3| , but recently it was shown that this claim is based 
on an erroneous analysis of the Lyapunov spectrum fll4| . 
However, for these two systems, no attempt of a theoret- 
ical understanding about the nature of an approximate 
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CPR has been carried out so far. In this paper, we ad- 
dress and attempt to clarify these issues. We find that 
for a dilute gas of hard spheres obeying the SLLOD dy- 
namics, where the masses and the radii of the spheres 
are not necessarily the same, and the total peculiar ki- 
netic energy is kept constant by applying the iso-kinctic 
Gaussian thermostat in between collisions, the CPR is 
violated at the most at 0(7 4 ), for small 7, in the ther- 
modynamic limit. [ jT5| | Our analysis is based on the key 
idea that the coefficient of friction representing the iso- 
kinetic Gaussian thermostat for a dilute gas of particles 
mutually interacting by means of a short-ranged poten- 
tial and obeying the SLLOD dynamics with a small shear 
rate, in the NESS, reaches a fixed value in the thermo- 
dynamic limit, with I/a/ZV fluctuations, where N is the 
number of particles. Eh]. 

The structure of the paper is the following: in Sec. II, 
we describe equations of motion for the SLLOD dynam- 
ics, define the Lyapunov exponents, and discuss the suffi- 
cient conditions for an exact CPR. In Sec. Ill, we demon- 
strate how the coefficient of friction representing the iso- 
kinetic Gaussian thermostat is expected to reach a fixed 
value at the thermodynamic limit, with l/>/~N fluctua- 
tions. In Sec. IV, we present the explicit calculations and 
discuss the status of an approximate CPR. To make the 
calculations In Sees. II-IV simple, we assume that each 
of the gas particles has a unit mass. Finally, we end 
this paper with discussions on possible generalizations, 
including the generalization to the case when the masses 
of the gas particles are arbitrary, in Sec. V. 



II. THE SLLOD EQUATIONS OF MOTION FROM 
A VIEWPOINT OF DYNAMICAL SYSTEMS 



Sill ( F » • Pi ~ JPixPiy) 
2-<i=l Pi 



a 



(2.2) 



The SLLOD equations of motion, without the dissipa- 
tive term — orpj, cannot be derived from a Hamiltonian 
(unlike the DOLLS equations of motion). 

We will use the equations of motion exclusively in 
terms of the particles' positions 17 and laboratory ve- 
locity Vj. This introduce s th e change of variable from 
to Vj = Pi +7?/iX. in Eq. (2.1), which can then be written 
as 



v 4 = P» 



(2.3) 



In the present context, the gas particles are hard 
spheres of arbitrary radii. This reduces the dynamics 
of the gas particles to an alternating sequence of flight 
segments and instantaneous binary collisions. During a 
flight, the dyna mics of the gas particles is therefore de- 
scribed by Eqs. (2.2 - 2.5 ) with Fj = 0. As far as the dy- 
namics in collisions is concerned, it is possible to derive 
the limiting behavior of the iso-kinetic Gaussian ther- 
mostat as Fj — ► 00 ]Tij ], but this leads to rather com- 
plicated collision rules. For the purpose of simplicity, 
in this paper, we choose to apply the iso-kinetic ther- 
mostat only in between collisions. Thus, at an instan- 
taneous collision between the i-th and the j-th sphere 
= 1, 2, .., N; i 7^ j), the post-collisional positions and 
laboratory momenta (+ subscripts) are related to their 
pre-collisional values (— subscripts) by 



{(v 4 _ - Vj_) ■n.jjni 



(2.4) 
(2.5) 



A. Equations of Motion 

The SLLOD equations of motion describe the dynam- 
ics of a collection of N particles constituting a fluid with 
a macroscopic velocity field u = 72/x (i.e., the gradient 
of the a;-component of the macroscopic fluid velocity u 
in the y-direction is 7). For simplicity, each gas parti- 
cle is assumed to have a unit mass. The specific form 
of the SLLOD equations of motion for the i-th. particle, 
in terms of its position r\ and peculiar momentum p^, is 
given by 



and 



r* = Pi + 72/i* > Pi 



JPry> 



api 



(2.1) 



where Fj is the force on the i-th particle due to the other 
particles in the system. The peculiar velocity of a parti- 
cle is defined as its velocity with respect to the velocity 
of the flow at its location and the peculiar momentum 
of a particle is the product of its mass and its peculiar 
velocity. The value of a, the coefficient of friction repre- 
senting the iso-kinetic Gaussian thermostat in Eq. (2.1), 
is chosen such that the total peculiar kinetic energy of 
the system, J2iPi/2, is a constant of motion, i.e., 



{0 



13 > 



(2.6) 



while the positions and the velocities of the rest of the 
spheres remain unchanged. Here, is the unit vector 
along the line joining the center of the i-th sphere to 
the j-th sphere at the instant of collision. Note that al- 
though in any particular collision, the peculiar kinetic en- 
ergy changes over a collision, these changes are random, 
both in magnitude and sign, due to the randomness of 
the collision parameters, and hence it is quite likely that 
the system would reach a steady state, where the average 
change of peculiar kinetic energy would be zero. 

To study the SLLOD dynamics as a dynamical sys- 
tem in three dimensions (the dimensionality does not af- 
fect our analysis), we form the 3iV-dimensional vectors 
R = (17, r 2 , ...,rjv), V = (vi,v 2 , Vjv) and Ny, whose 
l-th entry is given by N y = (c5 M - Sij) ny/v^ (I = 
1,2,..,N). Using these new variables, we write the 
SLLOD equations of motion in the fa, v^) coordinates 
of N hard spheres during a flight, Eq. ( |2.3| ), in a com- 
pact form 
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R 



V = a 7 CR - «V . 



(2.7) 



Here, C is a 37V x 3./V matrix with NxN entries, each of 
which is a 3x3 matrix. In terms of the entry index (I, m), 
in the xyz-basis, Q m = c<5; TO (l,m — 1,2, ..,N) and 



c = xy = 



1 





(2.8) 



At a collision between the i-th and the j-th sphere, the 
equations of motion are given by |ll[| 



R^ 



R 



V+ = V_ - 2 (V_ 



(2.9) 



In our ana lysis he reafter, except for Sec. Ill, we will use 
only Eqs. (2.7-2.9) to describe the dynamics. 

At this point, we introduce the following notations. A 
phase space point can be denoted as r = (R, V) . A 
linear transformation on phase space can be given as a 
6./V x 6-/V matrix. Any such matrix P can be split in 
terms of four 3N x 3N sub-blocks, for which we use the 
notation PW, pPl, P® and PW, such that 



p[l] p[2] 
p[3] p[4] 



(2.10) 



Each sub-block pW (i = 1, . . . , 4) itself can be divided in 
3x3 sub-blocks again, where each sub-block can be iden- 
tified by two indices I and m, I along the horizontal and 
m along the vertical direction (l,m = 1, 2, .., N). Such a 
sub-block is denoted by Pj^. 

B. Lyapunov Exponents for Hard-sphere Systems 

To calculate the Lyapunov exponents for hard-sphere 
systems, let us say that the system starts at time to 
at a phase-space location T(to) = (R(to), V(to)). Un- 
der time evolution, T(t) follows a trajectory in the 6N- 
dimensional phase space, which we call the "reference 
trajectory" . The set of N hard spheres would suffer 
a sequence of binary collisions on this trajectory. We 
also consider an infinitesimally displaced trajectory in 
the phase space, which starts at the same time to, but 
at r'(t ) = r(i ) + ST {t ). Under time evolution, T'(t) 
follows another trajectory, always staying infinitesimally 
close to the reference trajectory. This trajectory we call 
the "adjacent trajectory". We also assume that the set 
of N hard spheres on the reference and the adjacent tra- 
jectories suffer the same sequence of binary collisions. 
We denote the time evolution of the infinitesimal 6N- 
dimensional tangent vector 8T(t) over time (t — to) by 
the 6Nx6N matrix L(t — to), i.e., 



6T(t) = L(t - to) 5T(to) . 



(2.11) 



The Lyapunov exponents are the possible exponential 
growth rates in time of |l_r| for different directions of 



unit vectors f . We have to define the norm in an ap- 
propriate way. Making the time it takes for a sphere 
with a typical velocity vo to cross the distance of a typ- 
ical radius of a sphere ao our unit of time (i.e., ao/vo 
is set to 1) solves the problem that the components of 
r have different dimensions. For the inner product be- 
tween two tangent vectors 5T^ = (Sw- 1 ', SV^ 1 ') and 
6T {2) = (c5R (2) ,<5V (2) ) we use 



(5r (1) |<5r 



(2)\ _ 



5r 



(2) 



(2h 



The norm is now defined as |<£T| = -J (ST\ST) . 

The Lyapunov exponents are the logarithms of the 
eigenvalues of the matrix A, defined by 



A = lim 

t — >oo 



L(t-*o) 



l/[2(t-t )] 



(2.12) 



where, L(i— to) = [L(£— ^o)] T L(< — to). The corresponding 
directions for the exponential expansion and contraction 
of the phase space (R, V) are obtained from the eigen- 
vectors of L(f — to). 



The dynamics of ST(t) in Eq. (2.11), for a gas of hard 
spheres, can be decomposed into an alternating sequence 
of flights and instantaneous binary collisions. We denote 
the transformation of ST(t) over a flight segment between 
t and t + At by H(At) such that 

ST(t + At) = H(Af) 5T(t) with H(0) = I . (2.13) 



Explicitly, H(At) is obtained from 
ST(t) = T(t)ST(t) 



H(At) = exp T 



t+At 



dt! T(*') 



(2.14) 



(2.15) 



where the subscript T indicates time ordering. Notice 
that H(Ai) in a general system will depend on time t 
as well, but we have suppressed that in our notation. If 
we now denote the transformation of ST(t) over an in- 
stantaneous binary collision (say, between the i-th and 
the j-th sphere) by the matrix M^ , we can express the 
matrix L(i — to) in terms of the H and matrices in 
the following way: if the dynamics involves flight seg- 
ments separated by s instantaneous binary collisions at 
ti, ti . . .t s such that to < ti < t2 < ■ ■ ■ < t s < t, then 



L(t-to) = H<:A/jM, H(At._i).. 

...M njl H(At ). 



Here, AU = ij+i — U for i 
At«=t- t«. 



1,2,. 



(2.16) 
(s — 1) and 
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C. The Sufficient Conditions for an Exact CPR 



then also 



If the CPR is exactly satisfied for a dynamical system, 
the sum of the conjugate pairs of the Lyapunov exponents 
is some constant c, i.e., if Xi is a Lyapunov exponent of 
this system, then c—Xi is also a Lyapunov exponent. The 
proof of a possible conjugate pairing rule will follow from 
the properties of the matrix L(t— to). However, to under- 
stand the interplay between the properties of the matrix 
L(t — to) and an exact CPR in full generality below we 
first look at the property of L(t — to) that has been used 
in various cases to prove CPR. 

(a) If the matrix L(t — to) is symplectic, i.e., L(t — to) 
satisfies the symplectic condition 

[L(t-t )] T JL(t-t ) = J, 



Det[L(t- t ) - ft- 1] = 0, (2.20) 



for an eigenvalue L of the matrix L(t — to). That 
means that if L is an eigenvalue of L(t— to) then so is 
|U 2 L _1 . In that case, one finds from Eq. ( |2.12j ) that 
c =limt_ >00 (ln /i)/(t — to). If the system is ergodic, 
then this long time average for c can be equated 
to a NESS average. Notice that condition (a) is 
obtained as special case of condition (b), namely 
when fj, = 1. 



with J as the usual symplectic matrix, then 
L(t - to) J L(t - t ) = J. 



(2.17) 



Eq. (2.17) can be used to show that 

Det[L(t-t ) - LI] =0 
(with I the identity matrix) implies 



Det[L(f-f )-^l] =0. 

L 



(2.18) 



This means that if L is an eigenvalue of L(t — t ) 
then so is L _1 . Since the Lyapunov exponents are 
the logarithms of th e eigenvalues of the matrix A, 
defined in Eq. ( 2.12 ), it is easy to see that c = 0. 
All Hamiltonian systems fall in this class. 

(b) In the existing literature P-LT|, the concept of the 
symplectic condition defined above, has been gen- 
eralized to the so-called "^-symplectic condition" 
and applied to thermostatted systems where an iso- 
kinetic Gaussian thermostat keeps the total lab- 
oratory kinetic energy constant and the external 
force on the constituent particles of the system is 
dependent only on the positions of the particles. 
For these systems, in an appropriate reduced phase 
space characterized by all the non-zero Lyapunov 
exponents, the matrix L(t — to) satisfies this [x- 
symplectic condition, which means that there exists 
a time- dependent positive scalar quantity /i, such 
that [L(t - t )] T JL(t - t ) = mJ- This implies that 



L(i - 1 ) JL(t - 1 ) 



(2.19) 



which can be used to derive that if 
Det[L(t - to) - LI] = 



Returning momentarily to the SLLOD dynamics, we ob- 
serve that the formalism developed in Refs. [p|-|llfl fails 
here. The primary reason is associated with the fact 
that the total peculiar kinetic energy is held constant for 
the SLLOD dynamics, as opposed to the total laboratory 
kinetic energy in Refs. PHll]]- One however needs to in- 
terpret the statement regarding the connection between 
the violation of condition (b) and the non-exactness of 
the CPR with care. By virtue of the fact that condition 
(b) above is a sufficient condition for the CPR to hold ex- 
actly, the violation of an exact CPR cannot be guaranteed 
if condition (b) is not satisfied. 

Guided by this observation, the interplay between the 
properties of the matrix L(t — to) and an exact CPR for a 
dynamical system can be generalized further than what 
is presented in (b). If there exists any constant non- 
singular matrix K satisfying K 2 cx I and the following 
condition 



[L(t-t )] T KL(t-t ) =/iK 



(2.21) 



is sa tisfied with a time-dependent scalar quantity ^t, then 
Eq. ( 2.2C ) can be shown to hold for an eigenvalue L of the 
matrix L(t — to), implying that the CPR is exactly sat- 
isfied for such a dynamical system |]l9]| . In analogy with 



the nomenclature presented in (b), we call Eq. (2.21) a 
"generalized /i-symplectic condition" with matrix K. We 
emphasize that the necessary condition for an exact CPR 
to hold for a dynamical system is not known. 

In view of Eq. ( fOl| ), thus, one should look for such 
a matrix K to prove an exact CPR. Instead, we look at 
it from a different angle, namely that we would like to 
understand how the SLLOD dynamics of N hard spheres 
with an iso-kinetic Gaussian thermostat deviates from an 
exact CPR. 
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III. BEHAVIOR OF a IN THE 
THERMODYNAMIC LIMIT 

Our procedure to study the deviations from CPR be- 
gins with the following observation: in the thermody- 
namic limit, for the SLLOD dynamics with short-range 
inter-particle potentials at low density of spheres and at 
small 7, the behavior of a simplifies to a great extent. 
After some transient time the system reaches the NESS, 
and the coefficient of friction a fluctuates with l/^/N 
fluctuations around a fixed value ao |l6| . For not too 
large fluctuations, the distribution function for a can also 
be shown to be approximately Gaussian. Thus, to calcu- 
late the Lyapunov exponents for large iV at low density 
of spheres and at small 7, to which we con fine ourselves 
henceforth, a can be replaced by ao m Eq. (2/7). We will 
now briefly present the gist of the derivation in Ref . [jl6| , 
applied to hard spheres. 



For a hard sphere system, the force term in Eqs. (2.1) 
and (2.2) is zero during a flight. Thus, for a flight, we 
have 



■7 



' N 

E 



N 

^ PixPiy ■ 



(3.1) 



Introducing a second thermostat variable, 

-1 



JV 



Er 2 

i=l 



(3.2) 



a closed set of equation follows from Eq. (2.1) 



a = - let 
(3 = - 2a(3 . 



(3.3) 



These equations are valid during the flights, i.e., the in- 
tervals between collisions. 

We treat collisions by looking at their net effect, i.e., 
how the velocities and positions of the particles i and 
j involved in the collision, are changed from their pre- 
collisional values p^_ and Pj_ to their post-collisional 
values Pi+ and Pj+ . These are the only two velocities to 
change, and because Yli=iPi ^ s o f or der iV. the changes 



in a and (3 are, according to Eqs. (3^) and (3J5), of order 
TV -1 . The number of collisions in the whole system is an 
extensive quantity as well, so there are 0(N) of these 
0(1 /N) changes in a unit of time. The averages of the 
small changes are not zero, so there is a net effect of 0(1) 
per unit time to the time derivatives of a and /?, which 
we will denote by a respectively, b: 



n = 



2a 2 



$ = -2aj3 + b. 



(3.4) 



This set of equations has a fixed point (ao,/3o)j which is 
stable if a > 0, so the system reaches this fixed point 



after some time. On top of this dynamics, there are fluc- 
tuations. Assuming that the collisions are independent, 
the central limit theorem applies, and the fluctuations 
are 0(1/ \/N). For more detailed analysis, we refer to 
Rcf. Kq|, where the independence of the changes in a 
and P is linked to the assumption of molecular chaos. 

Finally, we note that to maintain a stationary total 
peculiar kinetic energy in a system with a constant ao 
thermostat, this constant has to be chosen differently for 
varying 7. For 7 = 0, i.e., in equilibrium, ao can be set 
to zero, and kinetic energy is determined by the initial 
conditions. Near equilibrium, i.e., in the linear response 
regime, the right hand side of Eq. fl3.l|) should scale as 7 2 , 



i.e. ao oc 7 . Obviously, there arc higher order correc- 
tions to this behavior which play a role for larger values 
of 7. If, for 7^0, the initial condition is such that the 
total peculiar kinetic energy is not equal to the stationary 
value, that value will be approached in time. 



IV. STATUS OF AN APPROXIMATE CPR IN 
THE THERMODYNAMIC LIMIT 

Based on the discussion in the last paragraph of Sec. II 
and using the results in Sec. Ill, we will explore the possi- 
bility of an approximate CPR for the the SLLOD dynam- 
ics of hard spheres in the thermodynamic limit, at small 
7 and at low density in this section. We will first obtain 
the desired results for a constant coefficient of friction ao 



in the equations of motion (2.7) and then we will discuss 
the validity of an approximate CPR when the system is 
under an iso-kinetic Gaussian thermostat. To this end, 
our starting aim is to study the generalized ^-symplectic 
properties of th e m atrix L(t — to) for the d ynam ics de- 
scribed by Eq. (2/7) during a flight and Eq. (2.9) during 
an instantaneous binary collision between the i-th and 
the j-th sphere. 

However, as the matrix L(t — to) is constructed from 
the H and the M matrices, we will have to study the 
generalized /i-symplectic properties of the H and the M 
matrices separately. 



A. Generalized /i-symplecticity Property of H(At) 

The matrix T describing the dynamics of 8T during 
flights is found from Eq. (2.7) to be 



T 



I 

ao7C — aol 



(4.1) 



From Eq. ( [4.l| ), it is straightforward to obtain 

H(At) = exp [TAt] , (4.2) 

where the 3x3 sub-blocks of H are given by hl|^(At) = 
hW( At) where 
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h«(At) 
h^(At) 

h( 3 >(At) 
h< 4 >(Ai) 



1 



7At- 

,— anAi 



7 [1 - cxp(- a At)] 



«0 



1 

a 

+X [aoAt(l 
a 

7 [1 

,-a Ai J | 



) - 2 + 2e 



e -«oA*l c 



and 



7 



At 



ao 



a Atl 



(4.3) 



Due to the complicated form of H(At), it is easier to 
study its /i-symplecticity properties in terms of the ma- 
trix T. This involves the task of finding a possible matrix 
K satisfying the condition 



T K + KT = /3K 



(4.4) 



such that K 2 oc I. If such a matrix K exists, then H(At) 
is generalized [i symplectic with that matrix, and 



jjL = exp 



t+At 



pdt' 



(4.5) 



Since K and T are constant ma trice s in the present con- 
text, (3 is also a constant. Eq. (4.4) can be treated as a 
simple eigenvalue equation to solve for the eigenvalue (3 
and the eigenvector K. We find that there exists a matrix 
G satisfying 



T T G + GT = -a G 



(4.6) 



and the 6Nx6N matrix G in terms of its 3 x 3 sub-blocks, 
is given by G^ = G^ = and Gj 2 ^ = — G^ = g, where 



1 

1 
1 



(4.7) 



We note that th ere may exist o ther forms of g such that 
G satisfies Eq. (4J3), but Eq. (4/7) is the simplest one 
that satisfies G 2 cx I and works for all 7. 

For the purpose of future use, we construct matrices To 
and H (At) by setting 7 = in Eq. (|tJ) and Eq. (|]|, 
without setting ao = (in reality, ao = when 7 = 0), 



i.e., 



To 
Ho (At) 



I 

-a l 

exp [T At] 



and 



(4.8) 
(4.9) 



More explicitly, the 3x3 sub-blocks of H (At) are given 

and h^^At) can be 



h ; (At) Si, T 



by (Ho)Pi(At) 
found by putting 7 = in Eq. (13) without putting 
ao = 0. The matrix Ho (At) is now not only generalized 
/i-symplectic with matrix G, but also /i-symplectic with 
J, i.e., 



[H (At)] T JH (At) = e-" oAt J. 



(4.10) 



The relevance of this observation will become clear in 
Sec. IV.C. 



We also note that Eqs. (4.1-4.7) hold for any constant 
coefficient of fri ction (not necessarily ao), which implies, 
using Eq. (2.21), that for a collisionless gas of point par- 
ticles obeying the SLLOD dynamics with a constant co- 
efficient of friction, the CPR is exact, as can be seen in 
simulations ]13[ ]. 



B. Generalized ^i-symplecticity Property of M;j 

Unlike the H matrices, the My matrices corresponding 
to a binary collision between the i-th and the j-th sphere 
do not follow from Eq. (2.9) directly. This is due to the 
fact that even though the sequence of binary collisions 
are the same on the reference and the adjacent trajecto- 
ries, the binary collision between the i-th and the j-th 
sphere on these two trajectories in the phase space are 
not simultaneous. One therefore needs the dynamics of 
the tangent vectors for the time interval, St, between the 
two collisions on the reference and the adjacent trajec- 
tories involving the i-th and the j-th sphere. To obtain 
an expression of My, we follow the formalism developed 
in Refs. |pr| , p2[ , which in turn, is based on the formalism 
presented in Refs. |@ and 

The dynamics of tangent vectors at a collision is de- 
rived in the Appendix. The result is that 



ST+ = My jr. 



with 



M, 



(I 



R 



I 



(4.11) 



(4.12) 



and R a symmetric matrix given by Eq. ( A18). This form 
of Mjj immediately implies that M is symplectic, but not 
generalized /i-symplectic with matrix G for [i = 1, i.e., 



but 



+ G. 



(4.13) 



(4.14) 



C. Generalized /i-symplecticity Property of L(t — to) 
and the Origin of an Approximate CPR 

From Sees. IV. A and IV. B above, we can finally see 
that for a collection of hard spheres obeying the SLLOD 
equations of motion with constant coefficient of friction 

(a) the H matrices are generalized /i-symplectic with 
matrix G, but not with matrix J (see Eq. (4.6)) 
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(b) the M matrices are symplectic but not g enera lized 
/ i-sym plectic with matrix G (see Eqs. ( 4.13 ) and 



jV 



Once the H a nd th e M matrices are combined together, 
following Eq. (2.16), the matrix 



L(f - f ) = H (At s )M ia j s H ( At s _ % ) . . . M 2lJ1 H(At ) (4.15) 

is seen to be generalized /i-symplectic neither with ma- 
trix G nor with matrix J. 

To study the degree of deviation fro m an exact CPR 
using the properties of L(t — t ) in Eq. ( 4.15 ) with con- 
stant coefficient of friction ao, we can use either K = G 
or K = J. While the former choice implies that one has 
to try to estimate the deviation from an exact CPR using 
the distribution of the unit vectors Ny's and the collision 
angles for different sets of binary collisions in the expres- 
sion of M, the latter choice means that one can make 
the estimate by using the property of the H matrices in 
Eq. (4.15). We choose the latter approach, because not 
only is it much easier to calculate the typical magnitude 
of a flight time of a sphere at low densities, but also, 
an estimate of the deviation from the exact CPR can be 
made at small 7, as a power series expansion in 7. How- 
ever, the smallness of 7, which has a dimension of inverse 
time, has to be defined in a proper manner. To do so, 
we notice that the density of the spheres n sets a time 
scale in the form of the mean flight time To of an individ- 
ual sphere, and in three dimensions To ~ hao/vy. Here, 
n = nciy is the dimensionless density and, as before, v 
and ao are typical velocity and radius of a sphere. Thus, 
the actual dimensionless small parameter corresponding 
to the shear rate is 7 = 7"r . 

A naive way to estimate the deviation from an exact 
CPR using the latter approach is to use the deviation of 
the H(At ) matrices from an exact /x-symplecticity (see 
Eq. ( 4.15| )). Such a deviation is characterized by the ma- 



trix D(At) = [H(At)] T JH(At) - e - Q ° A *J. The matrix 



D(At) can easily be calculated from Eq. ( [4.3|) . However, 
to estimate the order of the matrix elements of D(At), an 
order of estimate of the quantity At has to be obtained. 
To this end, we note that while To is th e me an flight time 
for an individual sphere, At in Eq. (4.15) denotes the 
mean time for a flight of N spheres. This implies that 
At ~ To/N, as on an average, there are N/2 different bi- 
nary collisions over a mean flight time To of an individual 
sphere. Thus, one would expect that in the thermody- 
namic limit, D(At) — ► D(0) = and o ne wo uld be led to 
conclude that the H matrices in Eq. (4.15) a re al l sym- 
plectic. This in turn would imply, from Eq. ( 4.15 ), that 
L(t — to) would be /i-symplectic and therefore a gas of 
hard spheres, obeying the SLLOD dynamics with a con- 
stant coefficient of friction a would satisfy an exact CPR 
in the thermodynamic limit. We demonstrate below that 
this simplification is not correct. 

The proper estimate of the deviation from an exact 
CPR has to be made by considering H(t ). To see why 
this is so, we rewrite the matrix H(At) as 



H(At) = e 



■ a At 



]JW(At) 



(4.16) 



with H^At) defined by 

(H l )M(At) = ^4<S M h( fe >(At) 

+ (l-M(4,i+<Me- QoA * /2 l]. (4-17) 

In effect, IT(At) describes the evolution of the infinites- 
imal deviation of the trajectory of the i-th sphere, while 
it has an almost trivial action on the infinitesimal devia- 
tion of the trajectory of the j-th sphere, j 7^ i. It is easy 
to see that H* (At) has some useful properties 

H*(A<i) H J (At 2 ) = H*(Ati + At a ) , 
[H i (Atj), H j (Atj)} = and 

[ H'(At), M icjc } = : .f i, / / and j, / i . (4.18) 



The properties of I T (At ) in Eq. ( 4.18| ) allow us to shuf- 
fle the terms in Eq. ( 4.15 ) such as to collect together as 
many H l (At)s with the same i as possible. The result is 
that to the right of any M, cJc figure an H lc (r?) and an 
W c (Tj a ), where t£ and t£ are the time of flights for the 
i c -th and the j c -th spheres before their mutual collision 
c. Consequently, 



TV 



L(t-to) = e— 



■ ao(t-to) 



n H '(*-^) 



n^^io^^). (4.i9) 



The product sign in Eq. (4.19) is to be expanded towards 



the left, i.e. Y\ c =i = A s ■ ■ ■ A\. Here, tj is the last time 
that i-th p article collided (or to if it didn't collide). From 
Eq. (4.19), it is now clear that the proper estimate for 
the deviation from an exact CPR has to be made by con- 
sidering the properties of H l (Tj), with n = 0(tq), and 
not from the properties of \~\ 1 (to/N). 

We also notice that if one uses the corr espon ding Ho 
matrices instead of the H matrices in Eq. (4.15) to con- 
struct an analogous matrix l_o(t — to), defined by 



L (t - t ) = H (At s ) M ls3s H (At s _!) . . 
...M hjl Ho(Ato), 



(4.20) 



then the matrix l_o(t — to) is /i-symplectic, because of 
Eqs. (4.1C) and (4.13). As a consequence, the logarithms 
of the eigenvalues of l_o(t — to), defined by 



Lo(t-to) = [Lo(t-to)] T L (t-to), 



(4.21) 



pair exactly. If we arrange the corresponding Lyapunov 
spectrum 
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A = lim 

t— too 



U>(t-to) 



l/2(t-t ) 



in the decreasing order of magnitude as X\ 



(4.22) 



(o) A (o) 



,A 6Ar , tnen + A 6Ar _ i+1 - -a . 



v(0) 



Th e ma trix E J (ri) can be easily calculat ed from 

Eqs. 04.2$ . With the aid of E q s - @) and ( EJli, we 
now express the matrix AH*(rj) defined in Eq. (4.23) as 

AH l fe) (^) = Ah J [ fe ](r 2 )^, m , where 



Motivated by this, our approach to study the devia- 
tion from an exact CPR for the matrix L(t — to) will be 
to take l_o(i — to) as the reference matrix. This is rea- 
sonable because, as we will show, L(t — to) and l_o(t — to) 
are very close for small 7, if t — to = r = 0(tq). To 
show this, we first write the matrix l_o(i — to) m the same 
form as Eq. ( |4.19| ), with H 1 (tj)s replaced by H (T t )s. We 
then relate, for t — t,Q = t = 0(to), the difference matrix 
AL(t— 1 ) = L(t — to) — \-o(t— to) to the difference between 
H«(ri) and Ufa). 

We now define new matrices AH 1 and E l , such that 

AW( n ) = H'fa) - Hj(Ti) 

E*(rO = [H^tOI^AH*^) (4.23) 

with 

(Ho)^(rO = ^[^h^ ) (r i ) 

+ (l-<5 ;!l )(4,i+4,4)e- Q0T ' /2 l] (4.24) 

[for the defin ition of hp' '* (r^), see the preceding paragraph 
of Eq. ( 4.1 0| ) ] . It is easily seen that E l (r^) has non-zero 
elements only for entries involving the i-th sphere. This, 
together with relations (4.18Q, implies that 



Ah l «(r,) 



7 r i 



7 [1 - exp(- a Ti)} 



00 



AU^(n) = ^[a T l (l + e' a ^) 



Ti ] c and 

n + — [l-e c 
a 



(4.26) 



The crucial point is that since ao <x 77, i.e., the dissipa- 
tion is quadratic in the shear field, AH l (ri) in Eq. ( 4. 26] ) 
can be expanded in powers of 7 to obtain 

ya T? 



Ah^(n 
Ah l(3) (ri) = jocQTi c and 
Ah^( n ) 



2! 
3! 



2! 



(4.27) 



[E l (r,), H j o (t,)] = fori^jand 
[ E 1 (t,), M icjc ] = if i c ^ i and j c ^ i . 



(4.25) 



at the mos t leadi ng order. Thereafter, having combined 
Eqs. (|4 23D and d4.26H4.27] ), all elements of E J (r t ) are 
seen to have a prefactor of order -y 3 . This makes E l (Tj) 
small compared to 1. Hence, to first ord er, AL(t) = 
L(r) — Lo(t) can be found from Eq. ( 4.19 ) by keeping 
only the terms linear in E l (Tj): 



JY 



AL(r) 



*° T LI H o^ + to-ii)\[ E E '( T + *° - *i) II M = H o «) 



-EIl M « H o «) H o (^)x [E^J + EM<)] II ^(^W^) 

c=l a=l b=c+l y 



(4.28) 



We now shuffle all the E Ic matrices to the right, by successively exchanging the order of E* c and the next 
M 6 H o ( T i b ) H o ( T j b ) t0 its ri g ht - These wiU commute very often, because most collisions in the sequence will involve 
other spheres than this i c -th one. If they do not commute, we write 

E fe M 6 H?(7-£)H£(7£) = M 6 H?«)H?(r*) [m 6 H£(t*)H*(t* 

From that point on, we work with "[M fc H b H^ 6 ]" 1 E i M fc H i> H^ b " , and shuffle that to the right. Repeating this, we end 
up with L (r) on the left side again, and find, symbolically, 



E^H^H^tM 



(4.29) 



S 

[LoMl^ALM =^{n(MH H )- 1 } [E^(tD + E^(tO] |jj MH Ho| 

c=l 

N 

£ {n( MH o H o) _1 } E J (r + to - tj) {[] MH Ho} . 



(4.30) 
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Since r is of the order of a flight time of a sphere, s is 
of the order of N, but a typical sphere will have suf- 
fered only a few [a few here is roughly 0(1)] collisions. 
Consider a typical term in the sum over c. The MHoHo 
terms in the products are from collisions that involve the 
i c -th or the j c -th sphere, or spheres that were involved 
in collisions with the i c -th or j c -th sphere. The number 
of MH H will thus also be of 0(1). 

The matrix E l ° has non-zero elements only for the en- 
tries associated with the i c -th sphere. The multiplication 
by MHqHo yields more non-zero entries associated with 
other spheres, whose number however is of 0(1). There- 
fore, in the whole sum of s = O(N) terms in Eq. ( |4~3C| ) 
will yield for typical entries a result of order one (and not 
of order N as the sum over s oc N suggests). Since all 
elements of the E m atric es are p roport ional to 77 2 as can 
be seen from Eqs. (^23|) , §L2§ and (p~27l) , we conclude 
for t — to = t — 0(to) that 



L(t) = Lo(t) [I + 7 3 B], 



(4.31) 



where the matrix B is of order 1 in 7 and order 1 in 
N. Note that B contains higher powers of 7 as well. 
Because it involves c and contributions from collisions 
between spheres, B is not proportional to the identity I, 
so we cannot regard it as a simple scalar factor (in which 
case the exact co njugat e pairing would be easy to obtain 
again). Equation J4.3l| ) implies 

AL(r) = L(r) - L (r) 

= 7 3 [B T L (r) + Lo(r)B] + 7 6 B T L (r)B. (4.32) 



We can now see, from Eqs. ( 4.31 ) and ( 4.32 ) that the dif- 
ferences between L(r) and Lo(t), and between L(r) and 
Lo (t) are small, as anticipated on page |^, by a relative or- 
der 7 3 . Therefore the logarithm of the eigenvalues of L(r) 
and L (t) also differ by a term of order 7 3 in an absolute 
sense. If we now divide the logarithms of these eigenval- 
ues by the time t, we see that the finite time (for time 
t) Lyapunov exponents, calculated from Lo and from L 

(which we denote as (r) and A;(t) respectively, for 
i = 1, 2 . . . 6N) differ by a term 0(7 3 /t) = 0(77 2 ). Us- 
ing the fact that that A^(r)s satisfy the conjugate pair- 
ing rule, Aj(r)s too, will satisfy it up to 0(77 2 ). From 
the condition that 



6JV 



(iJV 



Anr)=5>(T)=-3JVa , 



(4.33) 



i=i 



[which can be easily verified from Eqs. (4.3) and (4.12)], 
we see that on average over the pairs, the deviation from 
CPR is zero. 

We make one further observation at this stage. The 
Lyapunov exponents (even the finite time once) are in- 
variant under 7 — > —7, so in a power series expansion 
of the individual Lyapunov exponents in 7 J2lj , the odd 
powers vanish. We therefore conclude that the logarithm 



of the eigenvalues of L(r) and Lo(t) must differ by a term 
of order 7 4 , and hence conjugate pairing of Aj(r)s will be 
valid up to the correction of the form 77 s . 

To extend these ideas for large (t — to) [and finally for 
t — ► 00)], we proceed in the following way. Notice that the 
matrices L(t — to) and Lo(t — to) are positive definite and 
symmetric. This allows us to express the two matrices in 
the form of L (t — to) = ex p(Ao) and L(t — t ) = exp(A), 
where for large [t — to), both the eigenvalues of Ao and 
A must behave ~ (t — to). In these terms, the differ- 
ence between the Lyapunov exponents for L(t — to) & n d 
Lo(t — to) is related to A — Ao. Since the difference be- 
tween L(t — to) and l_o(t — to) has an explicit prefactor of 
7 3 , so will A — Ao. Using the symmetry argument that 
the Lyapunov exponents have to be even functions of 7, 
we obtain 



A; + ^6N-i+l 



ao + 0( 7 7 3 ). i = 1, . . . 6iV (4.34) 



To explici tly extend the formalism developed in Eqs. 
( 4.19 - 4.32 ) to large t — to and thereby obtain a rela- 
tion between A^s and A^s, we need to concatenate a 
lot of L(t)s. These do not commute with each other, 
nor do they commute with Bs in general. This prevents 
an explicit demonstration of how the deviation is built 
up. For the largest and the most negative Lyapunov ex- 
ponent, it has been possible to show that they pair to 
— ao plus corrections of 0(77 3 ) by means of a kinetic 
theory approach |2^,|23|], based on the independence of 
subsequent collisions of a sphere. One expects that in 
subsequent time-intervals of O(ro), the L(r) matrices are 
not qualitatively much different from each other. There- 
fore, w e exp ect that the coefficient of the 0(77 3 ) term 
in Eq. ( 4.34| ), to be of the same order as that for a flight 



time t = 0(to), i.e. of the order of B, which is O(l). 

We now return to the discussion of an approximate 
CPR for the SLLOD equations of motion with an iso- 
kinetic Gaussian thermostat in between collisions. As 
discussed in Sec. Ill, in the non-equilibrium steady state 
and in the thermodynamic limit, the coefficient of dy- 
namical friction representing the iso-kinetic Gaussian 
thermostat for the SLLOD dynamics fluctuates with 
1/yN fluctuations around the fixed value ao- We would 
therefore expect that the Lyapunov spectrum for the 
SLLOD dynamics with an iso-kinetic Gaussian ther- 
mostat are given by that of L(t — to) plus terms of 
0(1/ VN). Co nsequ ently, the approximate conjugate 
pairing of Eq. ( 4.34 ) can be extended to this thermo- 
stat once we neglect the 0(1/ \/N) terms in the sum of 
the corresponding A^ and Aejv-i+i- 



V. DISCUSSION 

In this paper, we started with a collection of N hard 
spheres, each with unit mass and arbitrary radius. Next, 
we argued how the coefficient of friction representing 







the iso- kinetic Gaussian thermostat in between collisions, 
fluctuates around a fixed value ao with 1 / \f~N fluctua- 
tions in the NESS. Using the properties of the transfor- 
mation matrices for the infinitesimal phase space clement 
ST, we then showed that the CPR is expected to be vi- 
olated, at the most, at 0(7 4 ) for constant coefficient of 
friction a>o. The source of the CPR violation is basi- 
cally the ao7CR term in Eq. (JjTij), as it is clearly seen 
that CPR would have been exact in the absence of this 
term. Finally, we extended that result to the case when 
the coefficient of friction represents an iso-kinetic Gaus- 
sian thermostat. Note that to obtain the deviations from 
CPR for the constant ao thermostat, the number of par- 
ticles did not matter, so for that case the result holds for 
any number of particles, but with an iso-kinetic Gaussian 
thermostat between collisions, we need to take N to infin- 
ity to make the connection to the constant ao thermostat. 

In addition, as mentioned in the introduction of this 
paper, the condition that the mass of each sphere be 
unity, is not necessary. If we assume that the mass of the 
i-th sphere is mj, then one can obtain the same equa- 
tions of motion, Eq. (£0-2.2) and (2.7-2.9) in terms of 

v' = 



the primed quantities defined by v\ 



«,r' 2 . 



and N-j, defined as 



N-i = 




in; 



n. 



(/ = 1, 2, .., N). It is then straightforward, albeit lengthy 
and laborious, to see that our entire analysis goes through 
in te rms of these primed variables, once we use U' in 
Eq. flAjj ) with 




(5.1) 
(5.2) 



However, we must note that even though the analysis pre- 
sented in Sees. IV is not affected when the masses and 
the radii of the spheres are not necessarily the same, one 
should not allow extreme variations in masses and the 
radii of the spheres. For large variation of the masses, 
the system may phase separate into phases in which the 
mean flight times are different, which will invalidate the 
use of the Boltzmann equation in Ref. [|l6j (the Boltz- 
mann equation has been used to show the approach of 
a towards the constant value ao in Ref. [@). Further- 
more, the use of a typical mean free time in Sec. IV. B 
will not be possible for large variation in the radii of the 
spheres. Our analysis in this paper therefore, holds for 
somewhat limited variation of the masses and the radii 
of the spheres. 

Another possible generalization of our analysis can be 
carried out for the case when the gas particles interact 
with each other by means of a short-ranged, repulsive, 
inter-particle potential (an attractive potential may cause 
bound states) for a constant multiplier thermostat. The 



dynamics of the particles can be again decomposed into 
flights and "collisions" at low density. While the transfor- 
mation of the infinitesimal phase space volume element 
over flig hts will again be determined by the H matrices as 
in Eqs. (4.2-4.2) (and thereby have the same properties 



as above), the corresponding My matrices will not neces- 
sarily have a_similar form and properties as presented in 

Eqs. ((yjynf). For a given average value of the peculiar 

kinetic energy, the constant multiplier thermostat has to 
be chosen carefully (see the last paragraph of Sec. III). 
Having denoted the constant coefficient of friction again 
by ao, where ao oc 7 2 , is easy to see that the resulting 
matrix IVl.y can be decomposed into a sum of two matri- 



ces, M\f and Al\%, where M y 



(0) 



is exactly /x-symplectic 

but AM ^ 

is the culprit for the violations of CPR in the first place) 
had been absent, then My would have been /i-symplectic 



|24fl and AM tJ would have been zero, but with Eq. ( |2.7D in 
its present form, it is easily seen that AM tJ is oc aj ~ 7 3 
and the proportionality constant depends on the ratio of 
the time the two colliding particles spend in contact with 
each other to the time of a flight. This ratio is very small 
at low density of the gas and therefore the elements of 



AM jj- are small in comparison with ! 



,(o) 



Since the H ma- 



trices have a similar property, one can repeat the analysis 
of this paper for such inter-particle short-range repulsive 

,(0) 

U 



potentials using Ho and MJ- as the reference matrices 
and arrive at the same conclusion as Eq. ( 4.34| ). 

Moreover, on the basis of the fact that at low densities 
and short-ranged repulsive interparticle potentials (that 
do not give rise to bound states), the constant multiplier 
thermostat and an iso-kinetic Gaussian thermostat are 
equivalent in the thermodynamic limit fl6f , one would ex- 
pect that under such conditions, the Lyapunov exponent 
spectrum for the Gaussian thermostat would be practi- 
cally the same as that of the corresponding constant mul- 
tiplier thermostat, where the constant multiplier is once 
again proportional to 7 2 . Therefore, in view of the discus- 
sion in the previous paragraph, we can also expect 0(f) A ) 
deviations from the CPR for a Gaussian thermostatted 
SLLOD equations of motion in the thermodynamic limit, 
where the particles interact with each other by means of 
a short-ranged repulsive potential. 

Finally, we note here that verifications of our theory by 
means of computer simulation remains a challenge. The 
work for simulations is in progress at present, although 
our preliminary experience suggests that to retrieve the 
77 3 scaling of the deviation of the sum of the pairs of 
Lyapunov exponents from — ao is not an easy task. 
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Q(r_) 



R 

V - 2 (V_ 



N,-,)N,: 



Next we notice that 



dQ_ 



ST*_ - T+ Si 



(A2) 



(A3) 



where, 8T*_ is the infinitesimal phase space separation 
between the two trajectories at A and C is 



ST* 



ST_ 



T_5t 



(A4) 



APPENDIX: 

In this appendix, we will derive the transformation 
of tangent vectors at a collision. The essential element is 
that the collision does not happen at the same time on 
reference and on the adjacent trajectory. To understand 
the origin of this time lag between the binary collisions 
between the z-th and the j-th sphere on the reference 
and the adjacent trajectories, in Fig. [j], we have depicted 
an exaggerated schematic picture of the collisions taking 
place in the 6./V-dimensional phase space on the refer- 
ence and the adjacent trajectories at points A and C re- 
spectively. The points B and D show the corresponding 
positions of the adjacent and the reference points respec- 
tively when the binary collisions at A and C are taking 
place. Thus, the pre-collisional separ atio n between the 
reference and the adjacent points is AB = <5r_, while 
the post-collisional separation is DC = ST + . Using that 
r.;| and \rj+Srj+VjST — r,— Sr^ — Vj5r| both have to 



equal a,i 



':) > 



the time lag St between the two collisions 



at A and C, therefore, can be easily expressed as 
(5rj_ - 5rj_) ■ rijj 



8t = — 



(Vj-_ - Vi_) 



n, 



<5R ■ N 4 
V_ • N„ 



(Al) 



To obtain the expression of ST + , w e first express the 
transformation of (R, V) in Eq. (2.9) over a binary col- 
lision between the i-th and the j-th sphere in a matrix 
form 



J 



<9Q 



-2V_ 



I 

dR y 



N, 



Adjacent 
trajectory 




Reference trajectory 
FIG. 1. A schematic diagram of the collision dynamics on 
the reference and the adjacent trajectories. 



and r± describes the equations of motion, Eq. (2.7), right 
after(before) the collision at A, i.e., 



(A5) 



R± 






.v±_ 




cto7CR± — aoV± 



Having calculated the quantity 



6>Q 

ar7 



from Eq. (fO), 



; <9R_ 



I - 2N, ; N, , 



(A6) 



where each entry of th e m atrix on the r.h.s. of Eq. ( |A6| ) is a 3Nx3N matrix, the expressions for 5R + and <5V + can 
be obtained from Eq. (A3) as 



5R+ = SR - 2 (5R_ • N y ) N y , 



and 



SV+ = A • SOL + (I — 2N ii N ii )5V_ + (I - 2N ii N ii )[a 7CR_ - a V_](5r - [a 7CR+ - cv V + ]5r, 



(A7) 



(A8) 



where A 



-2V_ 



OR 



and SR*_ = (5R- + V_ St). At this point, we use Eq. (|2.9|) and obtain 
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(I - 2N ii N tf ) [a 0l CR- - a V_ ] St - [a o7 CR+ - a Q V+ }St = ~ 2a o7 (Ny ■ CR ) Ny <5r . (A9) 
Following Appendix B of , the term A • SR*_ can be expressed as 

A • <5R*_ = -2 [(V_ ■ «Ni,-)N« + (V_ ■ N^)5Ny] , (A10) 



r 



where 

SN tj = -j= (0, 0, . ., Snij, ... - Sna, . . .0) , (All) 

satisfying ny • dny = 0. This orthogonality condition 

between and Stiij also implies that Ny • <5Ny = 0. 
To obtain an expression for 5rijj , we need to take a look 
at Figs. H and Figure || describes, in the laboratory 
frame, the binary collision process between the i-th and 
the j-th sphere on the reference and adjacent trajecto- 
ries; the thick-lined spheres are on the reference trajec- 
tory whereas the thin-lined spheres are on the adjacent 
trajectory. Figure || describes the same binary collision 
process in the reference frame of the i-th sphere (with 
center C). In Fig. |^, the thick- lined j-th sphere (with 
center D) on the left depicts the collision situation on 
the reference trajectory and the thin-lined j-th sphere 
(with center E) on the left depicts the collision situation 
on the adjacent trajectory. 



(ai+ ajXrijj + 5n,j) 



8r + p St 




FIG. 2. Collision between the i-th and the j-th sphere. 
Thick-lined spheres are on the reference trajectory whereas 
the thin-lined spheres are on the adjacent trajectory. 



Clearly, in Fig. ^, the infinitesimal vector Dill is given by 

{v j --v i -)5r (A12) 



Sr* 



5ya 



and since the lengths of both the lines CD and CE are 
a, + a,j {fli and cij are the radii of the i-th and the j-th 
sphere respectively), we have 



1 



5r* 



(A13) 



(a.+ aj)(n 




5r - 8r,_ + ( p - p* ) 5x t^\— 



("j +a ;)n.; 

FIG. 3. Same collisions as in Fig. El in the reference frame 
of the i-th sphere. 



Let us define a 3N x 3N matrix U composed of TV x TV 
blocks of 3 x 3 matrices, such that, in terms of the block 
indices the only non-zero entries of U are 



U« - -Uij = -U 3 -i = U# = -I , (A14) 



where I is the 3x3 unit matrix. One can now write, using 
Eqs. (|A^[A14|) and Eq. (|Al]) that 



SNij = r 1 

\/2(ai + a,j) 



[ U • <5R + U • f _ St] 



From Eqs. (A10) and Eq. (A15), we finally have 
v/2 



A-SR* 



V2(ai + a.j) 



(U-V-)Nj, 
V_ • N« 



SR 



(U- V_)Ni 



(I - 2X,,N,, i W • SR 



N y (U ■ V_) - (N r V_)U + V -' U . V ~ N^ 

V _ ■ IN o n 



(A15) 



SR 



(A16) 
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where W is a 3Nx3N sy mmet ric matrix. Finally, using 
Eqs. dAib , (0-© and (|A16|) , the expression for My 
can be obtained as 



= (I - 2N tf N y ) 



I 
R I 



where 



R =W -2a 7 



Njj • CR 

V_ • N, :i 



N- N- ■ 



(A17) 



(A18) 
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